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Abstract. In Part I of this paper we described an equilibrium model of a jet in the gravitational field corresponding 
to the rigid-rotation region of the galactic disk. We used linear stability analysis to find the waveguide-resonance 
instability of internal gravity waves due to the superreflection of these waves from the jet boundary. In this part of 
the paper, we perform nonlinear numerical 2D and 3D simulations of the development of this instability. We show 
that the shocks produced by this instability in the ambient medium of the jet are localized inside a cone with a 
large opening angle and are capable of producing features that are morphologically similar to those observed in 
galaxies with active nuclei (NGC 5252 for example). 



1. Introduction 



In Paper I ( Afanasiev et al., 2007[ ) we used linear stabil- 



ity analysis to investigate the possibility of the formation 
of regular structures resulting from the development of 
hydrodynamic instabilities in conical jets including those 
with large opening angles. At the same time, the very ex- 
istence of such weakly collimated jet outflows from active 
galactic nuclei raises certain doubts — observed forma- 
tions are widely believed to be radiation cones and not 
mass outflows (see the "Introduction" section in Paper I) . 
From purely hydrodynamic viewpoint it is not quite clear 
why in such similar accretion-jet objects protostellar sys- 
tems develop highly collimated jets, whereas active galac- 
tic nuclei produce wide-angle cone outflows. Note that 
it is extremely difficult to determine from observations 
whether we are dealing with a highly collimated or large- 
opening-angle conical outflow. This is due to the fact that 
jets from active galactic nuclei are almost without excep- 
tion observed through high-intensity shocks produced by 
the intrusion of the outflow matter into the unperturbed 
ambient medium ( [Falcke et al., 19961 |Nagar et~al., 1999) . 
Hence a significant number of shock-excited ions is present 
along the line of sight in any case. 

In this paper we analyze a situation that is alternative 
to that addressed in Paper I. In our case, the jet outflow 
from the nucleus is collimated, but the nonlinear stage of 
the development of instability in the jet results in the lo- 
calization of the resulting shocks inside a wide cone. We 
believe that this effect can explain the formation of struc- 



tures observed in the vicinity of Seyfert nuclei (Z-shaped 
"arms", "arcs", and "arches" — see Paper I). Note that 
Hardee (1982), Falcke et al. (1996) and Lobanov et al. 
(2006) already pointed out the possibility of the forma- 
tion of such structures in the vicinity of active nuclei due 
to the interaction of radio-jet matter with the ambient 
medium in the boundary layer. 

Section [2] describes the technique of numerical simula- 
tions. Section [3] discusses the results of these simulations 
and proves that perturbations that build up in the jet as 
a result of waveguide-resonance instability produce in the 
ambient medium a system of intensely radiating nonlin- 
ear waves. Section|4]summarizes the main conclusions and 
gives the final comments. 



2. Technique of numerical simulations 

The stationary model that we use in this part of the work 
is quite similar to that described in Paper I. The only 
exceptions are the numerical parameter values that we 
give at the end of this section. 



2.1. Basic Equations 

In our numerical simulations of the dynamics of perturba- 
tions we use the following set of hydrodynamics equations 
in divergent form written in the spherical coordinate sys- 
tem (r, 9, ip): 
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dt 
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djpU) , , ™. dp /|V|2 
^^+dw(pC/V) = --+p(^— -— j , 



(2) 



(5) 



^+div(.pI^V).-|+py^ctg^, (3) 

If +div [V(i? + p)] = 

a?- 
where 

or r sm 6/ at* r sm H o(p 

p is the density; p, the pressure; V = {U, W, V), the ve- 

locity vector; E — p—- h the total energy, and 

2 (7-1) 

p 

the internal energy per unit mass. 



£ = 



(7-l)p' 



2.2. Numerical Scheme and Boundary Conditions 

We numerically integrate equation set (Il])-([5]) using the 
TVD-E scheme ( |Ryu et al., 1993 ) adapted to the spheri- 
cal coordinate system with a variable radial step. 

To analyze the dynamics of axisymmetric perturba- 
tions, we use a two-dimensional (2D) scheme (r, 9) with 
integration domain (ri„ < r < r^x and < < 7r/2) 
containing a total of N,. x Ng cells. To analyze the dy- 
namics of nonaxisymmetric perturbations, we use a three- 
dimensional (3D) scheme (r, 9, ip) with integration domain 
{rin ^ r < rex', < < 7r/2; < (fi < 27r) containing a 
total of N,. X Ng X N^p cells. 

It follows from linear analysis (Paper I) that the per- 
turbation wavelength increases linearly with distance, be- 
cause kr = const. Accordingly, in our numerical simula- 
tions we use a variable step in radial coordinate Ar^+i = 
e^^Ar.i, where Ax — (Inr-ea; — lnri„)/iVr. With Ar so de- 
fined, the number of cells per wavelength remains constant 
along the radius. We set constant integration steps in the 
9 and if directions, i.e., A6' — 7r/2iVe, A(^ = in/N^. 

To avoid distorsions in the (r, 9) plane due to the dif- 
ferences between the scheme velocities of the propagation 
of perturbations in the r and 9 directions, Ar^ and r^A^ 
cells should have equal lengths in these directions. I.e., 
Aa; = M and Nr = 2Ne (Inrgj; - lnri„)/7r. 

For unstable modes to build up in numerical simula- 
tions, a finite-width transition layer is required between 
the jet matter and the ambient medium (see Paper I). We 
set the extent of the transition layer by setting the size 
of the cell in the ^-direction. In this case, the thickness 
of the transition layer decreases with increasing Ng. The 
equilibrium distributions in the transition layer have the 
following form: 

C/.(r) = C/,(r)/2, 



Ps{r) 



Rpj 



Fig. 1. Radial dependences of relative perturbations (the 
solid bold, solid, dashed, and dashed-and-dotted lines cor- 
respond to p, U, p, and W, respectively) for 9 — 9j at 
different time instants: t = 0.5 (top) and t = 1 (bottom). 
Large amplitudes correspond to the nj — harmonic of 
the waveguide-resonance mode and lower amplitudes, to 
the Kelvin-Helmholtz mode. 

Ps{r) ^Pj{r) ^Pa{r) , 

where subscripts "j" , "a" , and "s" refer to the jet, ambient 
medium, and transition layer, respectively, and Uj and Us 
are the radial velocities of gas flow in the equilibrium state. 

We set the initial perturbation of the 9 component in 
the following form: 



sin(fcr Inr + rrnp) 




where A^, is the initial perturbation amplitude (A^ <^ 
Uj); kr, the dimensionless wavenumber, and to, the num- 
ber of azimuth mode. We set the scale factor S equal to 
0.2. 

We use the following boundary condition. 

— In the symmetry plane of the system (9 = 7r/2): sym- 
metric boundary conditions /(7r/2 — 0) = /(7r/2 -|- 0) 
for E, p, U , and V, and antisymmetric boundary con- 
ditions /(7r/2 - 0) = -/(7r/2 + 0) for W. 

— At the symmetry axis of the system (9 = 0): 2D scheme 
— symmetric boundary conditions /(— 0) = /(+0) for 
E, p, and U and antisymmetric boundary conditions 
/(-O) = -/(+0) for W; 3D-scheme — f{r,-0,ip) = 
fir,+0,ip + n). 

— At <y9 = and ip ^ 2tt: periodic boundary conditions 
- /(-OH /(2^ - 0), /(2^ 0) = /(+0). 

— At the inner (r = ri„) and outer (r = rex) radial 
boundaries: 



/(r,„-0) = /„(n„-0)+/(r,„+A) 



S/(ri„ -0) 



foir-e. + O) + fir,^ - X) 



Bf (nn + A) 
Bf{re^ + 0) 



, f{re.+0) 



f}2r2 + C/2(r)/4 ' 



Bf{rex - A) ' 

where Bf(r) is the amplitude function (envelope 
of perturbations) and A, the perturbation wave- 
length. Subscript "0" indicates the equilibrium values. 
According to the results of linear analysis (Paper I), 
the amplitude function at the initial time instant is 
Bf{r) — r^f . At subsequent time instants Bf can be 
computed by approximating the minima and maxima 
of perturbations in the grid cells. 



2.3. Parameter Values and ttie Transition to 
Dimensionless Variables 

The dimensionless parameter r — tdyn/trad that describes 
the intensity of radiative cooling is determined by the 
ratio of two time scales of the problem: the dynamic 



Afanasiev et al.: Formation of lonization-Cone Structures.. II. 



3 



Fig. 2. Time dependences of the logarithms of density per- 
turbation amphtudes: the curves of families 1 and 2 corre- 
spond to 6* = 9j and 6 ~ 1.56j, respectively. The middle, 
lower, and upper curves of each family on the upper panel 
correspond to r{Nr/2), r(37Vr/4), and r{Nr/4), respec- 
tively. The solid curves in the lower panel were computed 
with a 256 x 128 grid and the dashed curves, on a 128 x 64 
grid. 



Fig. 3. Radial dependences of relative perturbations (p 
(the solid bold line); U (the thin line); p (the dashed line), 
and W (the dashed- and-dotted line)) at different time in- 
stants for 9 — Oj/2: t ^ 2 (top), t = 3.2 (middle), and 
t — 7.1 (bottom). 

Fig. 4. Same as Fig. [3l but for 9 ^9j. 



time scale tdyn and radiative cooling time scale trad- Here 
tdyn = f sm-9j/cj is the time scale of the propagation of 
perturbations from the boundary to the symmetry axis of 
the jet (9 = 0), and t^ad = p/[(l - \)CKP^e^l'^\ is the 
time during which the energy of gas decreases by a factor 
of e due to radiative cooling. At r ^ 1 the effect of ra- 
diative cooling has a negligible effect on the dynamics of 
perturbations (Paper I). At t > 1 radiative cooling plays 
an important part in the evolution of unstable modes. 

We do not take into account the processes connected 
to viscosity and heat conductivity because the basic es- 
timates show that the ratio of the dynamical time to 
the characteristic time of such processes has the order 
Kn/M <C 1, where Kn is the Knudsen number and M is 
the Mach number. 

We performed simulations for 7 = 5/3, t = 2 or r = 4, 
jet half-opening angles 9^ — 5° or 9j = 10° and the jet-to- 
ambient medium pressure ratio of i? = 5 or i? = 10. 

In accordance with the results of Paper I, we set 
the Mach number of the outflow equal to M ~ 0.69 or 
M ~ 0.74. For the initial perturbations, we set kr = 16 
and A„ = 10~^ or A^, = 10~*. We find no qualitative 
differences between the results of the two series and there- 
fore below we describe the simulations with A^, = 10~^. 
The results of the series of the numerical simulations have 
shown that harmonics with kr — 16 grow most quickly, it 
is only these harmonics that we show in the figures. 

Note that we had to set such a relatively large jet open- 
ing angle because of the limited potential of the computers 
employed. For perturbations to be processed correctly, the 
half-opening angle must contain at least 10 cells of the 
computation domain, whereas we could not increase N0 
appreciably because of the limited RAM of the computers 
employed. 

Dimensionless boundaries of the computational do- 
main were determined by the fin — 1 and ff,x — val- 
ues. In accordance with Paper I, we perform our anal- 
ysis for the radii corresponding to the rigid-rotation re- 
gion of the galactic disk with the angular velocity ~ 
100 — 300 km s~^kpc^^ — const (which corresponds to the 
region inside the bulge). 

The outer boundary in dimensionless units is equal to 
r^x — (0.5 — 3) kpc and the inner boundary is times 
smaller, ri„ ~ (22 — 130) pc, depending on the specifics of 
the rotation curves of real galaxies. 

Concerning the transition to dimensionless velocities, 
let us point out that the balance of pressures in the jet and 



in the ambient medium (equation (4) in Paper I) allows 
equation (5) of Paper I to be rewritten as: 



1 5c/| ^(P^_^ a* 



2 dr 



(6) 



We further take into account that Vj oc r and the 
dependence of gravitational potential on radius, in accor- 
dance with Paper I, is 



to find from 



Ui = Qr 



Pa 
Pj 



flr\/W 



(7) 



(8) 



It is hence natural to normalize velocity to flr^x. In 
this case, the sound speeds in the jet and in the ambient 
medium are unambiguously determined by the dimension- 
less parameters M and R. We have for typical galaxies 
rirea; ~ (100 — 300) km s^^, and hence ([5]) implies that we 
must be dealing with a high- velocity jet (for the R adopted 
above, e.g., Uj = Silr) despite the subsonic nature of the 
flow in this jet (Uj < cj). This is consistent with the con- 
clusion that we made in Paper I that jet matter must be 
strongly heated by the radiation of the galactic nucleus. 

We determined the dimensionless time as t = tH.. We 
have = (0.1 — 2) x 10~"'^^rads~^ and this corresponds to 

t = (0.5 - 9) X 10^^ s = (0.2 - 2.9) x 10^ years. 

Finally, to pass to dimensionless density, we chose the 
value Pj(rin) 

number densities nj(ri„) = 35 cm"-^ 



6.7 • 10 ^^gcm ^, which corresponds to 



3. RESULTS OF NONLINEAR NUMERICAL 
SIMULATIONS 

3.1. Basic Comments 

When analyzing the evolution of perturbations, we could 
make use of higher spatial resolution for axisymmetric 
modes (2D) and hence reproduce higher-order effects. We 
therefore begin with the dynamics of pinch (m = 0) modes 
(Section l3.2p . and only then discuss three-dimensional he- 
lical (m — 1) modes (Section 13. 3p . 

Computation of the luminosity and surface density of 
the simulated structures is associated with the following 
problems. It is evident from the set of equations ([I])-([5]) 
that the real luminosity in individual emission lines can- 
not be computed because the set includes no equations of 
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Fig. 5. Same as Fig. |3l but for d = 1.2561, 




Fig. 6. Velocities of shock fronts at time i = 7.1: angles 
\.2h9j (the upper solid line) and 1.50j (the lower solid 
line). The dashed and dashed-and-dotted lines show the 
sound speed in the ambient medium (ca) and the jet ve- 
locity {Uj), respectively. 



photoionization balance. However, structures discussed in 
the following sections are observed inside the ionization 
cones both in Balmer and shock-excited forbidden lines. 
Our analysis of the simulation results yielded the relative 
emission-line luminosity of the features of the wave pat- 
tern by integrating function {p^A{T) - p^A(To))/p^A(ro) 
along the line of sight. In this definition one has to sub- 
tract the contribution of hot matter inside the jet from the 
total radiation. We therefore did not integrate the above 
function over the cells located inside the jet {9 < 9j), be- 



cause the ionized-gas temperature at these sites is so high 
that the gas cannot radiate in optical lines and most of its 
luminosity is contained at X-ray spectral domain. 

3.2. Results of 2D Simulation 

As we pointed out in Section [21 we performed a set of 
computations with various parameter values and differ- 
ent numbers of cells in the computational space: Nr — 
128, Ng = 64, 128, 256. Careful comparison of eigen- 
functions, obtained from the linear analysis, with results 
of an initial stage of numerical simulations, shows that 
two unstable modes developed in all the numerical sim- 
ulations performed: the surface Kelvin-Helmholtz (KHI) 
mode and the mai waveguide-resonance IGW mode of 
the u~ family. Unlike the latter, the Kelvin-Helmholtz 
mode proved to be very sensitive to the thickness of the 
transition layer between the jet and ambient medium, 
which is determined by the step of the grid cell size 
in the 0-direction. The stabilizing effect of the velocity 
shift is a well known fact (see, e.g., |Bodo et al., 1989| 
[Mustsevoi V.V. &: Khoperskov, 19911 ). Therefore the KHI 
mode contributed appreciably to the resulting perturba- 
tion pattern only on the Ng = 256 grid (see Fig. [T]), 
whereas its amplitude was negligible in all other compu- 
tations. 

Figure [2] illustrates the above-mentioned effect of the 
"thickness" of the transitional layer on the perturbation 
increment: increment grows with increasing spatial resolu- 
tion of the grid and asymptotically tends to the value de- 
termined from linear analysis of the discontinuous model. 

A characteristic feature of the evolution of perturba- 
tions in all the computational series performed is that ra- 
dial dependences of the amplitude envelopes of all quan- 
tities can be fitted fairly well by function r^f (see Section 
3 of Paper 1) . At the same time, at the stage of nonlinear 
saturation these dependences transform to the form r°'f , 
which describes equilibrium distributions (see Figs.[3H5]). 

During the nonlinear stage the development of 
waveguide-resonance instability of the main (rij — 0) 
IGWs mode produces shocks in the medium that sur- 
rounds the jet with approximately paraboloidal fronts in 
the real three-dimensional space. Figure [5] shows the ra- 
dial dependence of the r component of mass velocity at 
9 = 12.5°. The perturbations discussed here can be seen 
to be real shocks given the well-known fact that the dis- 
continuity surfaces in computations based on the TVD-E 
numerical scheme are always "smeared" over at least five 
spatial cells of the system ( |Ryu et al., 1993[ ). It is clear 
from Fig.[5]that the width of the shock in our case is equal 
to four to six cells, i.e., we indeed observe a discontinuity 
surface — a shock front — at this location. A comparison 
of the perturbation propagation velocity Vsh in the direc- 
tion normal to the shock surface with the sound speed Ca 
in the medium that surrounds the jet corroborates this 



^ I.e., a mode with no eigenfunction zero points between the 
jet boundary and symmetry axis — Uj = (see Paper I). 



Afanasiev et al.: Formation of lonization-Cone Structures.. II. 



5 



2D t= 2.00 2D t= 3.20 2D t= 7.10 




3D T= 5.08 3D T= 15.06 3D T= 28.45 




Fig. 7. Luminosity contours for the two- (figures at the top) and three-dimensional (figures at the bottom) simulations 
in the {x = r sm9, y — r cos 6) plane at different time instants. The slanted line shows the jet boundary {6 = dj). 
Darker shade indicates higher luminosity values. 



conclusion. As is evident from Fig. [SI the shock-front ve- 
locity is supersonic and our analysis shows that the Mach 
number of the shock front, Mgh = Vsh/ca, does not change 
with radius. We believe that this increase of the velocity of 
the shock fronts with increasing radius can be explained 
by the combined effect of two factors. First, the devel- 
opment of instability with the distance from the center. 
Second, the well-known effect of the acceleration of the 
shock as it propagates toward decreasing density, because 
the distribution of the gravitational potential ([7|) adopted 
in our model implies that p oc . 

Figure [7] demonstrates the evolution of the relative 
perturbations of luminosity normalized to its equilibrium 
value. The propagation of shocks from the jet boundary 
into the ambient medium is immediately apparent. 

Let us dwell on the following point, which is of great 
importance for our analysis. A characteristic feature of 
the systems with Newtonian-type gravitational potentials 
(protostars) is that perturbations born in the jet encom- 
pass the entire surrounding "atmosphere". This is due 



to impedance in the atmosphere being lower than in the 
jet: PaCa < PjCj. That is why we can assume that in a 
protostar-type accretion-jet system the jet and the ac- 
cretion disk interact through waves in the atmosphere 
dLevin et al., 1999D . In our model, we have a reverse situ- 
ation: paCa > PjCj, which prevents the penetration of per- 
turbations into the ambient medium. In this case, shocks 
propagate inside a limited cone with the half-opening an- 
gle of 9cone about the central jet. 

Figure [8] shows the dependence of the characteristic 
depth of such a penetration in latitude 6 {6 cone — dj ) as a 
function of time. As is evident from this figure, 9cone first 
increases rapidly with time and then reaches saturation. 
Note that this effect — localization of perturbations in 
the ambient medium inside a wide cone around the jet — 
does not depend on the mode type, because impedance 
is independent of the particular type of wave symmetry 
(symmetric or helical). Moreover, we argue that this effect 
is also independent of the physical mechanism of mode 
excitation. This means that perturbations of both vol- 



6 



Afanasiev et al.: Formation of lonization-Cono Structures.. II. 



ume resonance modes due to superreflection and Kelvin- 
Helmholtz surface modes decrease equally rapidly in the 
ambient medium with the distance from the jet boundary. 

As is evident from Fig. [U the opening angle of the 
cone overtaken by perturbations exhibits quasi-periodic 
damped oscillations. This behavior is due to the above- 
mentioned nonlinear superposition of two modes with dif- 
ferent physical buildup mechanisms. These modes have 
different frequencies and, correspondingly, beating shows 
up at the half-difference frequency superimposed on the 
carrier frequency equal to the half-sum of the mode fre- 
quencies. After reaching nonlinear saturation, the KHI be- 
gins to smoothly decay and so does beating, and subse- 
quent evolution is determined by the dynamics of the vol- 
ume resonance mode rij = 0. Quasi-periodic oscillations 
show up most conspicuously in the contours of radially- 
averaged perturbations drawn in the "latitude — time" 
plane (see Fig. [9|) . 



30 




10 



2 3 ^- b 6 7 

t 

Fig. 8. Time dependence of the opening angle of the con- 
ical domain overtaken by perturbations. The bold solid 
line shows the results of the simulations with a 256 x 128 
grid, = 10°, T = 2, and R = 10. The thin solid hue 
shows the results of the simulations with a 256 x 256 grid, 
9j = 5°, T = 4, and R = 10. The dashed line shows the 
results of the simulations with a 256 x 256 grid, Oj = 5°, 
r = 4, 7? = 5. 

Oscillations are due to yet another effect. Heating of 
gas in shocks changes the dispersion law. This results in 
the decrease of the flux of energy transferred to pertur- 
bations due to the development of instability, and the 
decrease of the perturbation amplitude. Subsequent ra- 
diative cooling causes the relaxation of gas to a close- 
to-equilibrium state, the dispersion law is restored, and 
the perturbation amplitude grows again. Both processes 
described above (beating and quasi-periodic variations of 
the dispersion law) are essentially nonlinear and difficult 
to decouple from each other. 




Fig. 9. Contours of radially-averaged distributions of tem- 
perature {T{t) - T(0))/r(0) (top) and the 6'-component 
of velocity W{t) (bottom). 



2D t= 7.1 R= 8.9 3D t= 28.4 R= 11.8 




-40 -20 20 40 -40 -20 20 40 
fix (deg) Sx (deg) 



Fig. 10. Cross sections of the volume-density cone at 
the constant-radius level. Dark shades correspond to 
maximum-brightness regions. Left-hand panel shows the 
axisymmetric mode in the 2D simulation (the luminosity 
cone obtained by rotating the computed domain about 
the 9 = axis) . The right-hand panel shows the first heli- 
cal mode according to the results of 3D simulations. The 
central circle corresponds to the jet boundary. The corre- 
sponding radius and time are indicated above each figure. 

3.3. Results of 3D Simulations 

In the case of 3D simulations we are forced to restrict our 
analysis to a 128 x 64 x 32 grid, which is rather coarse 
in terms of the 9 and if angles. Therefore the growth of 
perturbations due to the effect described in Sect ion [5T^ was 
much slower than in the case of 2D simulations and ceased 
at lower relative amplitudes. It is, however, safe to say that 
all the basic patterns of the evolution of perturbations are 
the same as in the two-dimensional case and therefore here 
we discuss only the aspects that differ. 

First, in the case of 3D simulations we analyzed the 
development of the first helical (m = 1) mode until the 
development of shocks. In the constant-r cross section the 
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Fig. 11. Sky-plane projection of emission- line cones for different tilt angles i between the line of sight and cone axis. 
The upper block of figures shows the results of 2D simulations (axisymmetric perturbations) . The lower block of figures 
shows the results of 3D simulations (the first helical mode of perturbations). The upper maps show the luminosity 
(surface density) maps and the lower maps correspond to radial-velocity fields in arbitrary scale. 



wave pattern of this mode has the form of a one-armed 
trailing (relative to the rotation of the phase pattern) spi- 
ral — see Fig. 1101 In three dimensions this pattern has 
the form of a slowly rotating (compared to fir) helical 
structure filling a wide cone around the central coUimated 
jet. 

We constructed a series of simulated wave-pattern 
surface brightness maps ("emission-line images", see 
Section 13. ip for different tilt angles i of the cone axis to 
the line of sight. We also constructed the velocity fields, 
i.e., the maps of luminosity- weighted line-of-sight gas ve- 
locities. Figure [Tl] shows examples of such maps separately 
for the symmetric and helical patterns. 

It is evident from the figures that the sky-plane pro- 
jections of the three-dimensional helical pattern depend 
on angle i, allowing both the Z-shaped and loop-shaped 
emission-line patterns to be obtained. The rotation of 



this phase feature on the radial-velocity map results in 
switching from "red" to "blue" shift to be observed as 
we pass from one arm of the Z-shaped pattern to another 
(Fig. 12). We pointed out in the "Introduction" to Paper I 
that such a velocity behavior is typical of some galaxies 
with Z-shapcd emission-line features. In our model switch 
of velocities is explained first and foremost by the rotation 
of the wave pattern inside the cone. 



3.4. Astrophysical Applications 

Thus our simulations show that instabilities of highly col- 
limated jets emerging from active galactic nuclei result 
in the propagation of shocks from the jet boundary into 
the ambient medium. These shock waves are located in 
the cone with a half-opening angle of 15° 40° outside 



8 



Afanasiev et al.: Formation of lonization-Cone Structures.. II. 



the jet. It is evident that in actual situations the ionized 
gas behind the shock should radiate intensively in optical 
lines (i?a, [OIII] etc.). A rather short-wavelength struc- 
ture forms — fcr 15 — 20 — and hence even a relatively 
small tilt of the symmetry axis of the cone considered 
with respect to the observer would produce a projection 
effect making the structure to appear as a continuous ra- 
diation cone with brighter regular features — "arcs" for 
axisymmetric wave modes or a Z-shaped pattern for the 
first helical mode — superimposed on it. 

In Fig. [TT]we sketchy reproduce the situation described 
based on the results of our simulations for three different 
tilt angles of the cone symmetry axis to the line of sight 
and show the distribution of luminosity integrated along 
the line of sight (thus the adopted model assumes optical 
transparency of gas). It goes without saying that in real 
galaxies one has to take into account dust absorption be- 
hind the shock fronts. Thus the images of ionization cones 
reported by Ferruit et al. (1998) and Quillen et al. (1999) 
exhibit extended dust lanes associated with regular struc- 
tures inside the cones. 

The proposed scenario of the formation of ionized 
cones makes it also possible to describe to a first ap- 
proximation the radial-velocity pattern observed in these 
cones. In the previous section we already explained ve- 
locity switching inside a Z-shaped emission-line structure 
due to a rotating helical wave. At the same time, the ve- 
locity fields of some Seyfert galaxies exhibit only "blue" 
Doppler shifts in one ionization cone and only "red" shifts 
in the diametrically opposite cone. A good example is pro- 
vided by the NGC 5252 galaxy whose velocity field was 
pubhshed by Morse et al. (1998) and Moiseev et al. (2000). 
Below we show that such a behavior can be explained by 
the development of an axisymmetric wave mode. 

Indeed if the tilt angle of the radio-jet axis to the line 
of sight a ~ 30 — 60°, then the cone nearest to the observer 
should be seen through shock fronts propagating virtually 
toward the observer and the farthest cone, against reced- 
ing shock fronts (see Fig. [T2|) . It is evident from Fig. [6] 
that the velocities of these fronts can be estimated as: 



(1.2-M.6)ca = (1.2 H- 1.6) 



(1.2 -M. 6) 



RM 



(9) 



RM 



We now substitute the numerical values quoted in l2.3l into 
([9]) to determine 

Vsh ^ (150- 630) fcms"\ 

which agrees well with the observed Doppler shifts (see 
references to Paper I). Moreover, the velocities of shock 
fronts increase linearly with radius (see Fig. [6]), implying 
that the Doppler shift of radial velocities should also in- 
crease with the distance from the nucleus. 

In Fig. [13] we compare the [OIII]-line image of the 
NGC 5252 galaxy (contour image in Fig. 1 from Paper I) 
and examples of arbitrarily scaled simulated images for 
the modes with m = and m = 1. 




Fig. 12. A diagram explaining the formation mechanism 
of the observed velocity field in the case of the develop- 
ment of an unstable axisymmetric mode. We adopted the 
shock fronts in the r — 6 plane from the results of our 2D 
simulations for dimensionless time instant t = 10. 

Note that our simulations were either axisymmetric or 
nonaxisymmetric. However, in real situations pinch per- 
turbations and helical modes in the jet develop simultane- 
ously. Nonlinear superposition of such modes can explain 
the morphology of the wave pattern and the velocity field 
observed in the NGC 5252 galaxy if we assume that the 
helical and axisymmetric (pinch) mode dominate in the 
inner and outer parts, respectively. 

4. DISCUSSION AND CONCLUSIONS 

— Results of our nonlinear numerical simulations qualita- 
tively confirm the conclusions based on linear analysis 
made in Paper I. 

— The development of the waveguide-resonance instabil- 
ity of internal gravity modes in a highly coUimated jet 
outfiow from the galactic nucleus may result (in the 
domain controlled by the gravity of the bulge) in the 
formation of a short- wavelength {kr ~ 15 — 20) pe- 
riodic system of shocks (Mach cones) in the matter 
surrounding the jet. 

— The shock structure mentioned above embraces a wide 
cone with a half-opening angle of 6'c — 15° — 40°. The 
opening angle of the cone is determined by particular 
parameters of the outfiow (first and foremost, by the 
ratio of the impedances of the jet and the ambient gas). 

— Given the intense radiative cooling of gas behind the 
shock fronts and projection effects, it can be assumed 
that the shock system discussed here must appear to 
the observer as a wide radiation cone with a wave pat- 
tern superimposed. In this case, helical modes deter- 
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mine the appearance of Z-shaped emission-line struc- 
tures observed in a number of Seyfert galaxies. 

Note that a number of authors already solved a sim- 
ilarly formulated problem (see 



e.g., 



Rossi et al., 2000 



|Lim &: Steffen, 2dOT| [Hardee, 2003[ ), however, they ana- 
lyzed an earlier stage of the interaction between the jet 
and the surrounding gaseous clouds. A "cocoon" does in- 
deed form behind the bow shock when a supersonic jet 
intrudes into the ambient medium. This effect was demon- 
strated many times in numerical simulations and one can 
see it on images made with the Hubble Space Telescope, 
e.g., in the central region of Mrk 3 ( [Capetti et al., 1999] ) 
or Mrk 78 (Whittle et al., 2004). However, on longer spa- 
tial scale lengths (r > 0.5 — 1 kpc), like in the case 
of NGC3516, Mrk 573, and NGC 5252 (see images in 
Paper I), no such shock is observed. This unambiguously 
proves that the initial injection of the jet (which formed 
deep in the inner region of the system in the accretion 
torus around the black hole) occurred rather long ago. We 
believe that in these objects the jet has already punched 
a channel in a more or less dense medium in the central 
region of the galactic disk, and the bow shock has bro- 
ken into the intergalactic medium (like this happens with 
galactic fountains). Since then, the gas of the cocoon must 
have relaxed (at least partially) and it must interact with 
the jet over the scale lengths mentioned above, leading to 
instability. On long scale lengths the jet propagates ballis- 
tically, however, its gas loses energy as a result of radiative 
cooling and there are no new mechanisms for the buildup 
of oscillations (because of the extremely low density of the 
ambient medium), implying that gas does not show itself 
in any observable way. We simulated exactly this situa- 
tion. Note the instability of this type may develop even in 
the absence of a sharp jet boundary. It is sufficient for the 
supersonic velocity difference to occur over less than one 
wavelength, and there will always be such perturbations. 

We understand that our model is only a first step to- 
ward the construction of the pattern of regular structures 
in ionization cones. When constructing a self-consistent 
model of such objects one must take into account at least 
the following factors: 

— Contributions of the flat (large-scale disk) and triaxial 
(bar) subsystems. 

— Rotation of the galactic nucleus. 

— Slow rotation of the matter of the conical jet about its 
symmetry axis. 

— Interaction between the jet and the ambient gas of the 
disk (in the case of a small angle between the jet and 
the disk plane. 

— The presence inside the shock of a "cocoon" produced 
by the supersonic intrusion of gas of the conical outflow 
into the dense gas of the disk. 

We believe the above factors to be of great importance 
for the dynamics of various concrete objects, however, 
they do not determine the essence of the phenomenon and 
therefore we ignored them in this paper. However, when 



discussing the formation of the observed structures in con- 
crete objects one has to take most of the above factors into 
account. 

In addition, to conclude our discussion, we point out 
that coUimated jet outflows from active galactic nuclei 
leave the bulge region after having acquired an essentially 
nonlinear, shock structure. Let us also add that near the 
bulge region (in the sense that the gravitational poten- 
tial is dominated by the spherically symmetric part) this 
structure must undergo significant changes due to rather 
sharp change of the law of the distribution of the gravi- 
tational potential and hence that of pressure and density 
inside the jet and in the ambient medium. 
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Fig. 13. Left-hand panel: an [OIII]-linc contour image of the NGC 5252 galaxy. Central panel: simulated velocity 
map for the pinch mode (m = 0) i = 60°. Right-hand panel: simulated luminosity map for the helical mode (m = 1) 
i = 60°. 



This figure "figl.jpg" is available in "jpg" format from: 



http://arXiv.org/ps/astro-ph/0702639vl 



This figure "fig2.jpg" is available in "jpg" format from: 



http://arXiv.org/ps/astro-ph/0702639vl 



This figure "fig3.jpg" is available in "jpg" format from: 



http://arXiv.org/ps/astro-ph/0702639vl 



This figure "fig4.jpg" is available in "jpg" format from: 



http://arXiv.org/ps/astro-ph/0702639vl 



